Fluctuation-stabilized marginal networks and anomalous entropic elasticity 
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We study the elastic properties of thermal networks of Hookean springs. In the purely mechanical 
limit, such systems are known to have vanishing rigidity when their connectivity falls below a critical, 
isostatic value. In this work we show that thermal networks exhibit a non-zero shear modulus G 
well below the isostatic point, and that this modulus exhibits an anomalous, sublinear dependence 
on temperature T. At the isostatic point, G increases as the square-root of T, while we find G oc T" 
below the isostatic point, where a ~ 0.8. We show that this anomalous T dependence is entropic in 
origin. 
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The stiffness of elastic networks depends on the me- 
chanical properties of their constituents as well as their 
connectivity which can be measured by the average co- 
ordination of nodes. Maxwell showed that a network 
of simple springs will only become rigid once the con- 
nectivity exceeds a critical, isostatic value at which the 
number of constraints just balances the number of in- 
ternal degrees of freedom 1]. This purely mechanical 
and athermal argument can be applied in order to un- 
derstand the rigidity of such diverse systems as amor- 
phous solids [2], jammed particle packings and emulsions 
[3, y and even some folded proteins |5!]. Interestingly, 
under-constrained systems that are mechanically floppy 
can become rigid when thermal effects are present. Per- 
haps the best known example of this is entropic elastic- 
ity of flexible polymers 6:]. Even a single, freely-jointed 
chain that is mechanically entirely floppy becomes elas- 
tic at finite temperature T: such chains resist extension 
with a spring constant that is proportional to T. At the 
level of networks of such chains, the macroscopic shear 
modulus also grows proportional to T 6] , although a de- 
tailed understanding of such networks also depends cru- 
cially on effects such as swelling in a good solvent [7|. 
Many systems, including network glasses |8l4l0l| and some 
biopolymer networks Ill4l3l| can be considered interme- 
diate between a purely mechanical regime well above the 
isostatic point, and a purely thermal or entropic regime 
below the isostatic point. However, very little is known 
about thermal effects of such systems near the isostatic 
point [li-[l^. 



Here we show that simple model networks, consisting 
of randomly diluted springs, can be stabilized by ther- 
mal fluctuations, even at low connectivity for which they 
would be floppy at zero temperature. Interestingly, we 
find that the linear shear modulus G exhibits anoma- 
lous temperature dependence both at and below the iso- 
static point. Specifically, we find that G oc T", where 
a < 1. This is surprising since one might have expected. 



in analogy with freely-jointed chains, that such networks 
would exhibit ordinary entropic elasticity (G oc T) be- 
low the isostatic point, as the mechanically floppy modes 
are excited thermally. Moreover, we flnd two distinct 
anomalous entropic elasticity regimes in the connectivity- 
temperature phase diagram, with the Maxwell isostatic 
point acting as a zero-temperature critical point (Fig.[T|). 
We perform Monte Carlo (MC) simulations on 2D 



spring networks that consist of A'^ 



nodes arranged 



on a triangular lattice that are connected by iVsp — zN/2 
springs, where z is the average connectivity (2 = 6 for the 
fully connected network). Periodic boundary conditions 
are used in all directions. To avoid the problem of net- 
work collapse [l9|, we consider two cases: one in which 
we keep the system area A fixed and treat the springs as 
'phantom' (i.e., we ignore steric interactions, and hence 
the springs are allowed to overlap), and one where we fix 
the system pressure P and prevent the springs from over- 
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FIG. 1. Schematic phase diagram of thermal networks in the 
T~z representation, where 'reduced T' is the ratio of the tem- 
perature to the spring energy and z is the connectivity, with 
critical connectivity Zc. In a manner reminiscent of quantum 
critical points d?!, IS-], we find a critical regime that broadens 
out for temperatures above the T — critical point. 
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FIG. 2. Shear modulus G of networks with A'^ = 3600 nodes connected by phantom (i.e., potentially overlapping) springs at 
constant area A = Ao, where Ao is the area of a fully connected network at T = (main plots) and self-avoiding springs at 
P — Q (insets), (a) G as a function of connectivity z at various temperatures T* = fesT/fespfg — 10~® (lowest data), 10^*, 10^'^, 
10~^, 10~^ and 1 (highest data). Solid line shows T — Q results, (b) G as a function of T* for z = 6 (triangles), z ~ Zc = 3.857 
(circles) and z = 3 (squares), (c) Scaling of the shear modulus using the form G — \t^z\^ ,"? (r*|Az|~ ), where Az = z — Zc, 
for r* < 10~^. In both systems the asymptotes and the constants are the same, with a = 1.4 and h = 2.8. 



lapping (self-avoiding springs). In both cases the system 
energy is given by 



G is then given by 



u = ^-fY.^k,-i.)\ 



(1) 



where £, 



'-u ~ \l -^ij + Vij i^ *^^ length of the spring that 
connects node i to node j (for all pairs of nodes con- 
nected by a spring), £o is the rest length and fcgp is the 
spring constant. In order to lower the connectivity of 
the system we set kg-p = for randomly chosen springs. 
For the phantom network this is identical to removing 
springs, while for the self-avoiding network this method 
has the advantage of computational efficiency over sim- 
ply removing the springs, since springs with fcgp = still 
contribute steric interactions and hence the nodes are es- 
sentially confined to a 'cell' by the surrounding springs. 

To find the critical (isostatic) point Zc, where the sys- 
tem becomes non-rigid at T = 0, we use a conjugate 
gradient algorithm to calculate the shear modulus G at 
zero-temperature for each z. For two dimensional net- 
works Zc = 4, although due to finite size effects this 
will be somewhat smaller for each A'' value studied [20|. 
We then increase T in steps for each system (in both 
the fixed-area and self-avoiding ensembles) and allow the 
systems to equilibrate, obtaining configurations under 
shear. We note that there is an additional critical point 
zp ~ 2.084 [21|, corresponding to the connectivity perco- 
lation threshold, below which there is no connected path 
through the network. The shear modulus vanishes at this 
point for all temperatures. 

In order to shear the systems, we use Lees-Edwards 
boundary conditions 22j. For springs that connect the 
nodes at the top of the box to those at the bottom, Xij 
becomes Xij + "fLy, where Ly is the height of the simu- 
lation box, and 7 is the shear strain. The shear modulus 



G 



Ad-f^' 
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where T is the free energy of the system. For T = 0, 
J" simply reduces to U. It is not possible to directly 
calculate T from MC simulations, and the method we use 
to calculate the linear shear modulus G (i.e., G at 7 = 0) 



23| 



is given in the supplementary material in Ref. 

In Fig.[2l^a) we show the dependence of the shear mod 
ulus G on z for a range of temperatures T for both phan- 
tom and self-avoiding networks. At low temperature, the 
shear modulus in both systems closely follows the zero- 
temperature behavior, decreasing as z is decreased from 
the fully connected network. Below the critical point Zc 
we find that the shear modulus deviates from the zero- 
temperature behavior, becoming non-zero for all finite 
temperatures. For z > Zc the shear modulus is largely 
insensitive to temperature, while for z < Zc, G depends 
strongly on T. For high temperatures, the shear modulus 
becomes increasingly insensitive to z and deviates from 
the zero-temperature behavior at increasingly high con- 
nectivities above z^, until eventually, when kgT ~ fcsp-^O' 
the thermal energy of the system is such that the network 
structure becomes unimportant. 

The different regimes of the dependence of G on T can 
be seen in Fig. E^b). At high connectivities the shear 
modulus remains almost constant as the temperature is 
increased, rising only as the thermal energy ksT ap- 
proaches the spring energy kspi^. As we approach the 
critical point, however, we find that the shear modu- 
lus, which will be at T = 0, shows an approximate 
T^'^ dependence at low temperatures. This anomalous 
temperature dependence is apparent over many orders 
of magnitude, and in fact corresponds to the system be- 
coming stiffer than expected at low T for ordinary en- 



tropic elasticity. As we increase the temperature fur- 
ther, in the self-avoiding spring networks we see this T^ " 
dependence give way to linear T dependence, while in 
the phantom spring networks we see a steeper T depen- 
dence, although it does not become linear. For z < Zc 
we find another anomalous regime with G oc T", where 
a ~ 0.8, at low temperature followed by linear T depen- 
dence at high temperatures in both phantom and self- 
avoiding networks. As we see these anomalous regimes 
in both types of network, we conclude that they are 
not driven by steric interactions nor the particular self- 
avoidance method used, but instead by the random net- 
work structure of these low z value systems. Consistent 
with this, if we remove bonds in such a way as to leave 
one-dimensional chains of springs (i.e., chains with z — 2) 
or honeycomb lattices (with 2 = 3) we find G oc T even 
at low temperatures, as one would expect for ordinary 
entropic elasticity 23]. 

The observed shear moduli can be well described by 
a scaling form analogous to that of the conductivity of 
a random resistor network [2J] that has also been suc- 
cessfully used to describe the shear moduli of athermal 
spring and fiber networks [2C 
scaling Ansatz is given by 



25|. For our system, this 



G= |Az|"^(T*|Az| 



(3) 



where a and b are constants and Az — z — z^. We 
find the best collapse of the data at low temperatures 
(T* = ksT/kspil < 10-5) for both the self-avoiding 
and phantom networks using the exponents a = 1.4 and 
b — 2.8, as shown in Fig. [21(c). This again demonstrates 
the three low temperature regimes, with almost constant 
G for z > Zc, G scaling with T°-* for z < Zc and G show- 
ing T^'"^ dependence as Az -^ 0. We note that, similar to 
our findings, a recent study of athermal fiber networks in 
two dimensions, with both filament stretching described 
by fcgp and bond bending described by stiffness k, found 
that the shear modulus scales with fcgp k^^^ at the criti- 
cal connectivity [20|. 

The non-zero shear modulus we find in these thermal 
networks below Zc arises from the entropic contribution 
to the free energy. The shear modulus can be broken 
down into its energetic and entropic parts as 



G 



_ 1 fd^U 

= Ge + Gs, 
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(4) 



where 5* is the entropy, and both Ge and Gs can be 
calculated during our simulation runs [23[ . We first show 
the ratio of the entropic contribution to G, Gs/G, against 
the connectivity z in Fig. [3l[a) for a system of phantom 
springs with constant area. At low temperature we see 
that Gs/G rises sharply as z approaches Zc before satu- 
rating to Gs/G ~ 1 below Zc, corresponding to a dom- 
inant entropic contribution. For z > Zc, the energetic 
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FIG. 3. (a) Ratio Gs/G of the entropic contribution Gs to 
the shear modulus G as a function of /S.z = z — Zc for a 
phantom network at T* = fcsT/fcsp^o- (b) Behavior oi fiGs = 
— {cPS/d'y'^)/A as a function of Az for the same systems as 
above. Results are for A = Ao and A'' = 3600. 



contribution Ge dominates, although the entropic con- 
tribution Gs becomes increasingly important at higher 
T. 

Figure [3l^ a) suggests that the behavior below the criti- 
cal point can be understood in terms of Gs alone. Thus, 
when considering the origins of the anomalous tempera- 
ture dependence of the shear modulus observed in Fig. [51 
it is instructive to look at the behavior of d^S/d^^ with 
temperature and connectivity. From Eq. ([J) it can be 
seen that for pure entropic elasticity (where G oc T) 
we should see d^S/d^^ ex r°. In Fig. [S^b) we show 
PGs = —{d'^S/d'y^)/A against connectivity for a range 
of temperatures in a system of phantom springs at con- 
stant area. As can clearly be seen, (3Gs diverges at the 
critical point at low temperatures before decreasing again 
below Zc, where it is clearly not constant with tempera- 
ture. In Fig. HKa) we show l3Gs against temperature for 
a range of connectivities. At the critical point we observe 
that (3Gs oc T~^/^ at low temperatures, before it begins 
to decrease with a smaller gradient at higher tempera- 
tures. Similarly, for z — 3 < Zc we find /3Gs oc T^"-^, 
where we again see an initial decrease in f5Gs, although 
here it reaches a constant value at higher temperatures 
(/3Gs oc T°). The high value of d^S/d-/^ at low temper- 
atures (as compared to the constant value at higher T) 
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FIG. 4. (a) PGs = -(9^5/^7^)/^ against reduced temper- 
ature T* for phantom spring networks at area A = Aq with 
z ~ Zc = 3.857 (red circles) and z — 3 (blue squares), (b) 
Shear modulus G against temperature T for phantom spring 
network with z = 3 and spring constant fcsp = cT, where c is 
a constant. Lines show linear fits. 



corresponds to the entropy decreasing more rapidly as the 
system is sheared. As noted previously, for honeycomb- 
like lattices and ideal chains we find pure entromc elastic- 
ity which would show pGs oc T" throughout [23|. Hence, 
we conclude that the dependence of the entropy on shear 
strain 7 at low temperatures arises from the disordered 
nature of the network, leading to the anomalous temper- 
ature dependence of the shear modulus. We note that 
we see qualitatively similar behavior of l3Gs with T at 
low temperature for self-avoiding networks, as one would 
expect from Fig. [2jb). 

A possible origin of this anomalous temperature behav- 
ior in sub-critical networks could be the internal stress a 
of the network, which in the phantom networks arises 
from the resistance to the tension the network is placed 
under in order to maintain its area. This tension can be 
shown to be proportional to the temperature [23|. As 
such, at low temperatures the shear modulus can be ex- 
pected, on dimensional grounds, to scale as G oc cr^fcg"", 
which would appear as G oc T'^kl~°' in our simula- 
tions. A similar anomalous dependence on stress was 
found in athermal networks with disordered molecular 
motors in Ref. [26| . Interestingly, if one takes the spring 



pected for freely-joined chains linking nodes, then pure 
entropic elasticity would be recovered, with G oc T and 
d^S/dj^ oc T°. However, if fcgp — cT, where c is a con- 
stant, it follows from Fig.|4lja) that the gradient of G with 
T will depend on the value of c. In Fig. Hljb) we show 
the shear modulus against temperature for systems below 
the critical point with fcgp = cT using a range of c values. 
Though all the systems show linear T dependence we do 
indeed see that as c is made smaller the shear modulus 
becomes smaller, until we set c to be smaller than ~ 10^ 
where the results converge. 

Our results demonstrate that there are two distinct 
regimes with anomalous temperature dependence of the 
shear modulus, as illustrated in Fig.[TJ In both cases, the 
dependence on T is sublinear. Thus, at low temperatures, 
this corresponds to an anomalously large effect of ther- 
mal fluctuations. The natural energy scale in our model 
is the spring energy fcgp^O' which can easily be much 
larger than the thermal energy, even at room tempera- 
ture. For protein biopolymers, for instance, it is expected 
that fcgp ~ E(P /to, where the diameter d is of the order of 
nanometers and the Young's modulus E can be as large 
as 1 GPa [271, l28| , and hence the spring energy for a seg- 
ment of length £q ~ lOOnm can easily be more than 10^ 
times larger than fc^T at room temperature [29|. Hence, 
for such systems, reduced temperatures T* in the range 
< 10^^ can be relevant and network-level thermal fluctu- 
ations can be much larger than expected based on naive 
entropic estimates. Importantly, such network-level fluc- 
tuations are almost always ignored in prior fiber network 
models and simulations, where either purely mechanical 
models have been used (see e.g. Refs. |2Cll25l. l30l432| ). or 
where hybrid mechanical models that include only single- 
filament fluctuations have been used 33J, |3J| . Finally, it 
is interesting to note that our phase diagram in Fig. [1] is 
reminiscent of other systems with zero-temperature crit- 
ical behavior, such as quantum-critical points 



mm 



As in such systems, in which the critical point is also 
governed by fluctuations other than thermal, we find a 
broad critical regime that fans out and extends for tem- 
peratures potentially far above T = 0. 
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constant fcgp to be proportional to T, as would be ex- 
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In this supplementary material, we present the method used to calculate the shear modulus, its 
energetic and entropic contributions and the pressure of model spring networks. We show that the 
pressure of a sub-critical network (z < Zc) held at constant area is proportional to the temperature. 
Finally, we show that the linear shear modulus of both a honeycomb lattice and a freely-jointed 
chain are linear in temperature. 



I. METHOD 

We perform Monte Carlo (MC) simulations on systems 
of two-dimensional spring networks. The networks con- 
sist oi N = Ti? nodes arranged on a triangular lattice 
and connected by iVgp — zN/2 springs, where z is the 
connectivity [z — Q for the fully connected network). We 
consider two cases: in one case we keep the system area 
A fixed (using NVT simulations, where T is the temper- 
ature and we note that in our case the volume V corre- 
sponds to the area A) and treat the springs as 'phantom' 
(i.e. we ignore steric interactions, and hence the springs 
are allowed to overlap). In the other we fix the system 
pressure P (using NPT simulations, and hence the sys- 
tem area is allowed to change) and prevent the springs 
from overlapping (i.e the springs self-avoiding). In both 
cases the spring energy is given by 
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where iij is the length of the spring that connects node i 
to node j (if i and j are connected by springs), (.q repre- 
sents the rest length and fcgp is the spring constant. The 
total system energy U is then given as 
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where the sum is over all springs m. In order to lower the 
connectivity z of a network, we set kgp = for randomly 
chosen springs. For phantom networks this is identical 
to removing springs, while for self-avoiding networks set- 
ting fcgp = has the advantage of computational effi- 
ciency over removing the springs. For fully connected 
triangular lattice networks the overlap test simply con- 
sists of calculating the area of the triangles formed by 
the springs and checking if it changes sign when a node 
is displaced. Springs with fcgp = will still contribute 
steric interactions which allows for this overlap test to be 
used even at low connectivity, and hence the nodes are 
essentially confined to a 'cell' by the surrounding springs. 



The systems are allowed to relax at temperature T in 
an MC equilibration run of > 2 x 10® MC cycles. For 
the self-avoiding networks we then perform a production 
run of > 4 X 10® MC cycles to calculate the average di- 
mensions of the simulation box. We finally fix the self- 
avoiding networks to their average area and equillibrate 
them using NVT simulations. 

In the phantom networks, where the area is fixed, we 
are effectively applying a pressure on the network which 
will vary with the temperature. We may calculate this 
pressure in a production run using the expression for the 
virial pressure given as [1] 
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where ks is the Boltzmann constant, /,„ — fcsp(^m 
is the force of spring m which has length £„, and 
denotes the ensemble average. 

In order to shear the systems, we make use of Lees- 
Edwards boundary conditions [2| in our equilibrated sys- 
tems. The energy of springs that connect the nodes at 
the top of the box to those at the bottom (via periodic 
boundary conditions) is modified to become 

k 



i^sp(^.,) = ^(y0^^7+T^7+4-^") ' (4) 



where Ly is the height of the simulation box, and 7 is the 
shear strain. The shear modulus G is given by 

where T is the free energy of the system. For T = 0, 
J- = lA and as such G is simply given by the second 
derivative of the energy of the system divided by the area. 
As it is not possible to directly calculate F from MC 
simulations we instead calculate the shear stress a from 
which we can obtain G. We consider the difference in 
free energy A J^ between a sheared and unsheared system 
(7 — 0) given as 

AJ- = -fcsT[ln(Z^) - In(Zo)] (6) 
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where Z is the partition function of the system, /3 = 
l/{kBT), and the sum is over all states k. The energy 
of the sheared system is U-y and Uq is the energy of the 
unsheared system. Taking the derivative of AJ-" with 
respect to 7 and dividing by the area A, we obtain the 
shear stress as 

_ 1 aj- _ 1 d/\T 
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giving the shear modulus as G = Ge + Gs- The system 
energy U is calculated as the ensemble average energy 
during the simulation run (U), which is defined as 
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(7) 



which is the ensemble average of the derivative of total 
system energy with respect to 7. This is simply equal to 
the sum of the derivatives of all spring energies 
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dj 



(8) 



We note that this derivative is non-zero only for springs 
with the energy given by Eq. (j4]) (i.e. those crossing the 
periodic boundary conditions). The derivative of Eq. ^ 
is given by 



dEsp^r,, ^'P^y \\lylj+^'^v+lLyY - ioj {x^j + -/Ly) 



dj 



\/yfj + (Xrj +lLyf 



(9) 

We calculate the average value of the Eq. ([9]) during sim- 
ulation runs for a range of 7 values, from which we obtain 
G simply as the derivative of Eq. ^ 
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(10) 



To calculate the entropic and energetic contributions 
to the shear modulus, we first note that the free energy 
can be defined as 



T = U-TS, 



(11) 



where S is the entropy. We can hence write the shear 
stress as 



_]_dU_TdS 
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from which we define 



1 dU 
TdS 



(13) 



By taking the derivative of this with respect to 7 we 
obtain 
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and with Eq. ([H]), Eq. ^^ and Eq. (O we can show 
1 



os = ^P 
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(17) 



We use this together with Eq. p4|) to obtain the entropic 
contribution to the shear modulus. 



II. PRESSURE OF CONSTANT AREA SYSTEM 

We calculate the pressure of a network using the virial 
pressure equation given in Eq. (jSj for connectivity z = 3 
and show the pressure P as a function of temperature T 
in Fig. [H As can be seen, we find that these networks are 
under tension (P < 0), with a magnitude that increases 
linearly with the temperature. 




FIG. 1. Pressure P of a phantom spring network with z = 3 
at constant area A = Ao (where Ao is the rest area of a 
fully connected network at T = 0) as a function of reduced 
temperature T* — ksT /ksp(.%. Points are simulation data, 
line shows a linear fit. 



and 



III. 



SHEAR MODULUS OF ORDERED 

SYSTEMS 




FIG. 2. Example configurations of fioneycomb lattice net- 
work (top) and freely-jointed chain (bottom). Blue lines are 
springs, red points are nodes. 



configurations of which are shown in Fig. [2]) using the 
method given in Section HI Both these networks are sub- 
critical (z < Zc) and each node in the networks has the 
same connectivity {z — 2 for the freely-jointed chain, 
2 = 3 for the honeycomb lattice networks). The behavior 
of the shear modulus G of these systems with tempera- 
ture T, along with that of a randomly diluted triangular 
lattice network, is shown in Fig. |31 where it can be seen 
that G shows linear T dependence in both the honeycomb 
and freely jointed chain systems, but not in the random 
network. 
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FIG. 3. Shear modulus G of honeycomb lattice network (HC), 
freely-jointed chain (FJC) and randomly diluted triangular 
lattice network (RN) as a function of reduced temperature 
T* — ksT /ksplio with area A = 0.9j4o, where Ao is the rest 
area of a fully connected triangular network at T = 0. Points 
show simulation data, lines are to guide the eye. 



We calculate the linear shear modulus G of both freely- 
jointed chains and honeycomb lattice networks (example 
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